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Abstract 

THERMUS is a package of C++ classes and functions allowing statistical-thermal 
model analyses of particle production in relativistic heavy-ion collisions to be per- 
formed within the ROOT framework of analysis. Calculations are possible within 
three statistical ensembles; a grand-canonical treatment of the conserved charges 
B, S and Q, a fully canonical treatment of the conserved charges, and a mixed- 
canonical ensemble combining a canonical treatment of strangeness with a grand- 
canonical treatment of baryon number and electric charge. THERMUS allows for 
the assignment of decay chains and detector efficiencies specific to each particle 
yield, which enables sensible fitting of model parameters to experimental data. 
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Nature of problem: 

Statistical-thermal model analyses of heavy-ion collision data require the calculation 
of both primordial particle densities and contributions from resonance decay. A set 
of thermal parameters (the number depending on the particular model imposed) 
and a set of thermalised particles, with their decays specified, is required as input 
to these models. The output is then a complete set of primordial thermal quantities 
for each particle, together with the contributions to the final particle yields from 
resonance decay. 

In many applications of statistical-thermal models it is required to fit experimen- 
tal particle multiplicities or particle ratios. In such analyses, the input is a set of 
experimental yields and ratios, a set of particles comprising the assumed hadron res- 
onance gas formed in the collision and the constraints to be placed on the system. 
The thermal model parameters consistent with the specified constraints leading to 
the best-fit to the experimental data are then output. 

Solution method: 

THERMUS is a package designed for incorporation into the ROOT [2] framework, 
used extensively by the heavy-ion community. As such, it utilises a great deal of 
root's functionality in its operation. Three distinct statistical ensembles are in- 
cluded in THERMUS, while additional options to include quantum statistics, reso- 
nance width and excluded volume corrections are also available. 
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THERMUS provides a default particle list including all mesons (up to the -fC|(2045)) 
and baryons (up to the Q~) listed in the July 2002 Particle Physics Booklet [3]. For 
each typically unstable particle in this list, THERMUS includes a text-file listing 
its decays. With thermal parameters specified, THERMUS calculates primordial 
thermal densities either by performing numerical integrations or else, in the case 
of the Boltzmann approximation without resonance width in the grand-canonical 
ensemble, by evaluating Bessel functions. Particle decay chains are then used to 
evaluate experimental observables (i.e. particle yields following resonance decay). 
Additional detector efficiency factors allow fine-tuning of the model predictions to 
a specific detector arrangement. 



When parameters are required to be constrained, use is made of the 'Numerical 
Recipes in C [1] function which applies the Broyden globally convergent secant 
method of solving nonlinear systems of equations. THERMUS provides the means 
of imposing a large number of constraints on the chosen model (amongst others, 
THERMUS can fix the baryon-to-charge ratio of the system, the strangeness den- 
sity of the system and the primordial energy per hadron). 



Fits to experimental data are accomplished in THERMUS by using the ROOT 
TMinuit class. In its default operation, the standard function is minimised, 
yielding the set of best-fit thermal parameters. THERMUS allows the assignment 
of separate decay chains to each experimental input. In this way, the model is able 
to match the specific feed-down corrections of a particular data set. 



Running time: Depending on the analysis required, run-times vary from seconds 
(for the evaluation of particle multiplicities given a set of parameters) to several 
minutes (for fits to experimental data subject to constraints). 
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1 Introduction 



The statistical-thermal model has proved extremely successful p!|2|3] in de- 
scribing the hadron multiplicities observed in relativistic collisions of both 
heavy-ions and elementary particles. The methods used in calculating these 
yields have been extensively reviewed in recent years [^115] . The success of these 
models has led to the creation of several software codes PI7|8] that use ex- 
perimental particle yields as input and calculate the corresponding chemical 
freeze-out temperature (T) and baryon chemical potential (/is). In this pa- 
per we present THERMITS, a package of C++ classes and functions, which 
is based on the object-oriented ROOT framework [9J. All THERMITS C++ 
classes inherit from the ROOT base class TObject. This allows them to be fully 
integrated into the interactive ROOT environment, allowing all of the ROOT 
functionality in a statistical-thermal model analysis. Recent applications of 
THERMITS include P|lUfllfl2fl3fl4fl5fl6fl7fl8fl9] . An on-going effort to ex- 
tend the range of applications of THERMITS has led to several publications 
on fluctuations in statistical models [20|21|22] . 

The paper is structured in the following way. In Section 2 an overview is 
presented of the theoretical model on which THERMITS is based. Section 3 
outlines the structure and functionality of the THERMITS code, while Section 
4 explains the installation procedure. 



2 Overview of the Statistical-Thermal Model of Heavy-Ion Colli- 
sions 

2. 1 Choice of Ensemble 

Within the statistical-thermal model there is a freedom regarding the ensemble 
with which to treat the quantum numbers B (baryon number), S (strangeness) 
and Q (charge), which are conserved in strong interactions. The introduction of 
chemical potentials for each of these quantum numbers (i.e. a grand-canonical 
description) allows fluctuations about conserved averages. This is a reasonable 
approximation only when the number of particles carrying the quantum num- 
ber concerned is large. In applications of the thermal model to high-energy 
elementary collisions, such as pp, pp and e^e~ collisions [23p^ . a canonical 
treatment of each of the quantum numbers is required. Within such a canon- 
ical description, quantum numbers are conserved exactly. In small systems or 
at low temperatures (more specifically, low VT^ values), a canonical treat- 
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ment leads to a suppression of hadrons carrying non-zero quantum numbers, 
since these particles have to be created in pairs. In heavy-ion collisions, the 
large number of baryons and charged particles generally allows baryon number 
and charge to be treated grand-canonically. However, at the low temperatures 
of the GSI SIS, the resulting low production of strange particles requires a 
canonical treatment of strangeness |25j. This is the so-called mixed-canonical 
approach. 

In order to calculate the thermal properties of a system, one starts with an 
evaluation of its partition function. The form of the partition function obvi- 
ously depends on the choice of ensemble. In the following sections, we consider 
the three ensembles most widely used in applications of the statistical-thermal 
model. 



2.1.1 The Grand-Canonical Ensemble 

This ensemble is the most widely used in applications to heavy-ion colli- 
sions |^26f27f28f29f3Uf31f32f33f34f35] . Within this ensemble, conservation laws 
for energy and quantum or particle numbers are enforced on average through 
the temperature and chemical potentials. 

In the case of a mult i- component hadron gas of volume V and temperature 
T, the logarithm of the total partition function is given by. 



lnZ«^(T,y,{M)= E 7f^/^'pln(l±e-^(^-^-))"\ (1) 



where gi and /i, are, respectively, the degeneracy and chemical potential of 
hadron species i, (3 = 1/T, while Ei = \Jp^ + , where mj is the particle 
mass. The plus sign refers to fermions and the minus sign to bosons. 

Since in relativistic heavy-ion collisions it is not individual particle numbers 
that are conserved, but rather the quantum numbers S, S and Q, the chemical 
potential for particle species i is given by, 

Hi = BifiB + Sifis + Qif^Q, (2) 

where Bi, Si and Qi are the baryon number, strangeness and charge, respec- 
tively, of hadron species i, and /i^, f^s and /ig are the corresponding chemical 
potentials for these conserved quantum numbers. 



species I 
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Once the partition function is known, the particle multiphcities, entropy and 
pressure are obtained by differentiation: 



5^^ = A(TlnZ«^), (4) 
P- = T^^. (5) 
Furthermore, the energy is given by. 



E^- = T^^-^^+ E f^^^r. (6) 

species i 

Using the prescription for the particle multiplicity. 



where we have introduced the z^. Similar expressions exist for the energy, en- 
tropy and pressure. 



In practice, the Boltzmann approximation (i.e. retaining just the k = 1 term 
in Equation ([7])) is reasonable for all particles except the pions. In this ap- 
proximation. 



lnZ«^(T,r,{/i,})= E ^,jd'pe-^^^^-^^^= E z]e^^^, (8) 

species i \^'^) species i 

where z] is the single-particle partition function of hadron species i. Further- 
more, under this approximation, P = I^species i N^'~^T /V both for massive and 
massless particles, which is certainly not true for quantum statistics. 



Since the use of quantum statistics requires numerical integration (or evalu- 
ation of infinite sums), while Boltzmann statistics can be implemented ana- 
lytically, it is worthwhile to identify those regions in which quantum statis- 
tics deviate greatly from Boltzmann statistics. In most applications of the 
statistical-thermal model, only a small region of the n — T parameter space is 
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of interest. Using the freeze-out condition of constant E/N [36], the thermal 
parameters, and hence the percentage deviation from Boltzmann statistics, 
can be determined as a function of the colhsion energy y/s. From such an 
analysis it is evident that, for pions, quantum statistics must be implemented 
at all but the lowest energies (deviation at the level of 10%), while, for kaons, 
the deviation peaks at between 1 and 2%. For all other mesons, the deviation 
is below the 1% level. For baryons, the deviation is extremely small for all 
except the protons at small y/s. 

When quantum statistics are applied, restrictions have to be imposed on the 
chemical potentials so as to avoid Bose-Einstein condensation. The Bose- 
Einstein distribution function diverges if. 



Such Bose-Einstein condensation is avoided, provided that the chemical po- 
tentials of all bosons included in the resonance gas are less than their masses 
(i.e. iii < nii). 



2.1.2 The Canonical Ensemble 

Within this ensemble, quantum number conservation is exactly enforced. Con- 
sidering the fully canonical treatment of B, S and Q in the Boltzmann ap- 
proximation, as investigated in [37], the partition function for the system is 
given by. 



(9) 




X exp 2 



mesons j 



E 




(10) 



where. 



E 



1 Qi(Sj^S+Qj(t>Q) 



baryons j 
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Zq represents the contribution of those hadrons with no net charges, and the 
sums over mesons and baryons extend only over the particles (i.e. not the 
anti-particles) . 



Once the partition function is known, we can calculate all thermodynamic 
properties of the system. Using thermodynamic relations it follows that. 



d 

S^-^{TlnZB,s,Q), 



fill 



and. 



dlnZB,s,Q 
dV ' 



(12) 



Furthermore, the multiplicity of hadron species i within this ensemble, N^''^''^ , 
is calculated by multiplying the single-particle partition function for particle 
i, appearing in the canonical partition function, by a fictitious fugacity Aj, 
differentiating with respect to Aj, and then setting A^ to 1: 



B,S,Q ^ din ZB,S,Q{>^i 



Ai=l 



(13) 



Following these prescriptions. 



B,S,Q 



^B-Bi,S-Si,Q-Qi \ j^GC 



(14) 



S^''''^ = lnZs,s,Q+ E ( 

species i \ 

pB,S,Q = ^ [ 
species i \ 

^B,S,Q ^ [ 

species i \ 



Zn- 



B—Bi,S—Si,Q—C 



species I \ B,S,Q 

z 



B-Bi,S-Si,Q-Qi \ pGC 



^GC 



Hi=0 



T 



B- Bi,S-Si,Q-Qi pjGC 



W=0 ' 



W=0 



(15) 
(16) 
(17) 
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One notices that, in the Boltzmann approximation, the particle and energy 
density and pressure of particle species i, within the canonical ensemble, differ 
from that in the grand-canonical formalism, with all chemical potentials set to 
zero, by a multiphcative factor {ZB-Bi,s-s^,Q-Q^/ Zb,s,q)- This correction fac- 
tor depends only on the thermal parameters of the system and the quantum 
numbers of the particle (i.e. the correction for the A+ and p are the same). 
The entropy is, however, slightly different; the total entropy cannot be split 
into the sum of contributions from separate particles. 



Now, 



lim 



'B-B,,S-S^,Q-Q^ 
^B,S,Q 



^B,t,B/T^S,t,s/T^Q,tiQ/T^ 



Thus, for large systems, the grand-canonical results for the particle number, 
entropy, pressure and energy are approached [37] . 



2.1.3 The Mixed-Canonical (Strangeness- Canonical) Ensemble 

Within this ensemble, the strangeness in the system is fixed exactly by its 
initial value of S, while the baryon and charge content are treated grand- 
canonically. For a Boltzmann hadron gas of strangeness S, 



27r 

— TT 

X exp 



5s e 



-iS<t>s 



E 



giV 



-hadrons i 

where the sum over hadrons includes both particles and anti-particles and 



(19) 



(20) 



Applying the same prescription for the evaluation of the particle multiplicities 
as discussed for the canonical ensemble, it follows that. 



Nf 



N 



GC 



(21) 



Furthermore, 
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S'' = lnZs+ E 

species i 

'Z 



s-s. 



E. 



GC 



GC 



species I 



species i 



S—S-i \ pGC 



MS=0 



(22) 

(23) 
(24) 



As in the case of the canonical ensemble, the strangeness-canonical results, 
in the Boltzmann approximation, differ from those of the grand-canonical 
ensemble, with fis = 0, by multiplicative correction factors which depend, in 
this case, only on the thermal parameters and the strangeness of the particle 
concerned. For large systems and high temperatures, these correction factors 
approach the grand-canonical fugacities, i.e.. 



lim 

V— >-oo V Zs 



(25) 



The expression for Zs can be reduced [38] to. 



+ 00 +00 



Zs = ZqX E ^|3m+2n-5|(a;i) /|n| (3:2) /|m| (a^s) 



(26) 



m=— 00 ?i=— 00 



where Zq is the contribution to the total partition function of the non-strange 
hadrons, while, 



Xi = 2 \ k^ik 



and. 



I k+i 
k—i 



with, 



(« = 1,2,3), 



(27) 



(^ = 1,2,3), 



(2J 



hadrons j with Sj =m 



MS=0 



V. 



(29) 



In |39fi0] it is suggested that two volume parameters be used within canoni- 
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cal ensembles; the fireball volume at freeze-out, V/, which provides the overall 
normalisation factor fixing the particle multiplicities from the correspond- 
ing densities, and the correlation volume, Vc, within which particles fulfill 
the requirement of local conservation of quantum numbers. In this way, by 
taking < V/, it is possible to boost the strangeness suppression. In fact, 
this was shown to be required to reproduce experimental heavy-ion collision 
data [39B0] . 



2.2 Additional Features 

2.2.1 Feeding from Unstable Particles 

Since the particle yields measured by the detectors in collision experiments in- 
clude feed-down from heavier hadrons and hadronic resonances, the primordial 
hadrons are allowed to decay to particles considered stable by the experiment 
before model predictions are compared with experimental data. For example, 
the total 7r+ yield is given by, 

iV^+= ^ iV,(P™)5r(z^7r+), (30) 

species i 

where Br{i — )■ 7r+) is the number of 7r"'"'s into which a single particle of species 
i decays. As shown in Figure [H approximately 70% of vr+'s originate from 
resonance decay at RHIC energies. Thus, a full treatment of resonances is es- 
sential in any statistical-thermal analysis. 



2.2.2 Resonances and the Inclusion of Resonance 
Width 

The inclusion of a mass cut-off in the measured resonance mass spectrum is 
motivated by the realisation that the time scale of a relativistic collision does 
not allow the heavier resonances to reach chemical equilibrium [H]. This as- 
sumes that inelastic collisions drive the system to chemical equilibrium. If the 
hadronisation process follows a statistical rule, then all resonances should, in 
principle, be included [31j. This is problematic, since data on the heavy res- 
onances is sketchy. The situation is saved by the finite energy density of the 
system, resulting in a chemical freeze-out temperature at RHIC of approxi- 
mately 160-170 MeV |12P3] . which strongly suppresses these heavy resonances 
and justifies their exclusion from the model. It is, however, important to check 
the sensitivity of the extracted thermal parameters to the chosen cut-off. 
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Fig. 1. The energy dependence of the decay contribution to the final tt~^ yield as pre- 
dicted by the statistical-thermal model. Calculations performed within the grand- 
canonical formalism, assuming the constant E/N freeze-out criterion |36) (weak 
decays excluded). 

The finite width of the resonances is especially important at the low tempera- 
tures of the SIS. Resonance widths are included in the thermal model by dis- 
tributing the resonance masses according to Breit-Wigner forms P5|24f27f25f41| 
This amounts to the following modification in the integration of the Boltzmann 
factor 125): 



J (fp exp 



d^p J ds exp 



f 



mV 



TT (s — m^)^ -|- m^r^ ' 



(31) 



where F is the width of the resonance concerned, with threshold limit ^^threshoid 
and mass m, and ^/s is integrated over the interval [m - 6m, m + 2F], where 
5m = mm[m - mthreshoid, 2F]. 



2.2.3 Deviations from Equilibrium Levels 

The statistical-thermal model applied to elementary e~^e~, pp and pp colli- 
sions [23112^ indicates the need for an additional parameter, 75 (first intro- 
duced as a purely phenomenological parameter [^5P6] ). to account for the 
observed deviation from chemical equilibrium in the strange sector. Since a 
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canonical ensemble was considered in these analyses, there is an additional 
strangeness suppression at work, on top of the canonical suppression. Although 
strangeness production is expected to be greatly increased in AA collisions, 
due to the larger interaction region and increased hadron rescattering, a num- 
ber of recent analyses [44ll47|48ll49ll50|51j have found such a factor necessary 
to accomplish a satisfactory description of data. 

Allowance for possibly incomplete strangeness equilibration is made by multi- 
plying the Boltzmann factors of each particle species in the partition function 
(or thermal distribution function fi{x,p)) by 75 ', where l^il is the number of 
valence strange quarks and anti-quarks in species i. The value 75 = 1 obvi- 
ously corresponds to complete strangeness equilibration. 

It has been suggested [22] that a similar parameter, 7^, should be included in 
thermal analyses to allow for deviations from equilibrium levels in the non- 
strange sector. Furthermore, as collider energies increase, so does the need for 
the inclusion of charmed particles in the statistical-thermal model, with their 
occupation of phase-space possibly governed by an additional parameter, 7c. 



2.2.4 Excluded Volume Corrections (Grand- Canonical Ensemble) 

At very high energies, the ideal gas assumption is inadequate. In fact, the total 
particle densities predicted by the thermal model, with parameters extracted 
from fits to experimental data, far exceed reasonable estimates and measure- 
ments based on yields and the system size inferred by pion interferometry [53] . 
It becomes necessary to take into account the Van der Waals-type excluded 
volume procedure [53|54f55j . At the same fixed T and fis, all thermodynamic 
functions of the hadron gas are smaller than in the ideal hadron gas, and 
strongly decrease with increasing excluded volume. 

Van der Waals-type corrections are included by making the following substi- 
tution for the volume V in the canonical (with respect to particle number) 
partition function, 

V^V- Yl (32) 

hadrons i 

where Ni is the number of hadron species i, and z/j = 4 (4/37rrf) is its proper 
volume, with its hard-sphere radius. This then leads to the following tran- 
scendental equation for the pressure of the gas in the grand-canonical ensemble 
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(assuming h particle species): 



P(T,/ii,..,/i,) = X:i^^'^'''(r,/I,), (33) 

2=1 

with, 

fli = fXi-UiP{T,jji,...,iJ,h). (34) 
The particle, entropy and energy densities are given by. 



ni T,/ii,...,/i?, = — — 35 

1 + Ej Ujnf''°-'{T, Hj) 



1 + Ej'^irif-'(T,/I,)' 



and, 



E,ef^"'(T,/I, 
l + Ei^^jnj.'^'"'(^,Aj)^ 



respectively. One sees that two suppression factors enter. The first suppression 
is due to the shift in chemical potential /Xj — )■ /Ij. In the Boltzmann approxima- 
tion, this leads to a suppression factor e'^^^^^ in all thermodynamic quantities. 
The second suppression is due to the [1 + Vjuf^^^T, fij)]'^ factor. 



In ratios of particle numbers, although the denominator correction cancels 
out, the shift in chemical potentials leads to a change in the case of quantum 
statistics. In the Boltzmann case, even these corrections cancel out, provided 
that the same proper volume parameter i/ is apphed to all species. 



3 The Structure of THERMUS 



3. 1 Introduction 



The three distinct ensemble choices outlined in Sections I2.1.H|2?T731 are im- 
plemented in THERMUS. As input to the various thermal model formalisms 
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one needs first a set of particles to be considered tfiermalised. Wlien combined 
witli a set of thermal parameters, all primordial densities (i.e. number density 
as well as energy and entropy density and pressure) are calculable. Once the 
particle decays are known, sensible comparisons can be made with experimen- 
tally measured yields. 



In THERMITS, the following units are used for the parameters: 



Parameter Unit 

Temperature (T) GeV 

Chemical Potential (/i) GeV 

Radius fm 



Quantities frequently output by THERMUS are in the following units: 



Quantity Unit 
Number Densities (n) fm~^ 

Energy Density (e) GeV.fm"^ 
Entropy Density (s) fm~^ 

Pressure (P) GeV.fm"^ 
Volume {V) fm^ 



In the subsections to follow, we explain the basic structure and functionality 
of THERMUS (shown diagrammatically in Figure [2]) by introducing the ma- 
jor THERMUS classes in a bottom-up approach. We begin with a look at the 
TTMParticle objectH] 



^ It is a requirement that all ROOT classnames begin with a 'T'. THERMUS 
classnames begin with 'TTM' for easy identification. 
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Particle 

TTMParticle 




Parameter Set 

TTMParameterSet. 



Partiole piopeities 
Partiole decays 



1 



Thermal Particle 

TTMThermalPartlcle 



Functions to calculate piimordial 
particle-, energy- and entropy 
densities and pressure 



Particle Set 

TTMParticleSet 



Thermal Model 

TTMThermalModel 



■ Primordial properties and feed-down 
contributions of all constituents calculable 

■ Total fireball properties now known 



plus experimental data 



Thermal Fit 

TTHThermaUit 



Best -fit thermal parameters can 
now be determined 



Fig. 2. The basic structure of THERMUS (only the most fundamental base classes 
are shown). 

3.2 The TTMParticle Class 



The properties of a particle applicable to the statistical-thermal model are 
grouped in the basic TTMParticle object: 

********* LISTING FOR PARTICLE Delta(1600)0 ********* 

ID = 32114 
Deg. = 4 
STAT = 1 

Mass = 1.6 GeV 
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Width 



0.35 GeV 
= 1.07454 GeV 



Threshold 



Hard sphere radius = 



B = 1 



S = 
Q = 

Charm = 
Beauty = 
Top = 



ISI = 



ICI = 



UNSTABLE 



Decay Channels : 



Summary of Decays : 



Besides the particle name, 'Delta (1600)0' in this case, its Monte Carlo nu- 
merical ID is also stored. This provides a far more convenient means of ref- 
erencing the particle. The particle's decay status is also noted. In this case, 
the A (1600)° is considered unstable. Particle properties are input using the 
appropriate 'setters'. 



3.2.1 Inputting and Accessing Particle Decays 

The TTMParticle class allows also for the storage of a particle's decays. These 
can be entered from file. As an example, consider the decay file of the A(1600)°: 



11.67 


2112 


111 


5.83 


2212 


-211 


29.33 


2214 


-211 


3.67 


2114 


111 


22. 


1114 


211 


8.33 


2112 


113 


4.17 


2212 


-213 


15. 


12112 


111 


7.5 


12212 


-211 



Each line in the decay file corresponds to a decay channel. The first column lists 
the branching ratio of the channel, while the subsequent tab-separated integers 
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represent the Monte Carlo ID's of the daughters (each hne (channel) can 
contain any number of daughters). The decay channel list of a TTMParticle 
object is populated with TTMDecayChannel objects by the SetDecayChannels 
function, with the decay file the first argument (only that part of the output 
that differs from the previous listing of the particle information is shown) : 



root [ ] part->SetDecayChaimels("$THERMUS/particles/Delta\(1600\)0_decay.txt") 
root [ ] part->List() 

********* LISTING FOR PARTICLE Delta(1600)0 ********* 



UNSTABLE 



Decay Channels : 



BRatio: 


0, 


.1167 


Daughters : 


2112 


111 


BRatio: 


0, 


.0583 


Daughters : 


2212 


-211 


BRatio: 


0, 


.2933 


Daughters : 


2214 


-211 


BRatio: 


0, 


.0367 


Daughters : 


2114 


111 


BRatio: 


0, 


.22 


Daughters : 


1114 


211 


BRatio: 


0, 


.0833 


Daughters : 


2112 


113 


BRatio: 


0, 


.0417 


Daughters : 


2212 


-213 


BRatio: 


0, 


.15 


Daughters : 


12112 


111 


BRatio: 


0, 


.075 


Daughters : 


12212 


-211 



Summary of Decays : 

2112 20% 

111 30.34% 

2212 10% 

-211 42.66% 

2214 29.33% 

2114 3.67% 

1114 22% 

211 22% 

113 8.33% 

-213 4.17% 

12112 15% 

12212 7.5% 
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In many cases, the branching ratios of unstable hadrons do not sum to 100%. 
This can, however, be enforced by scahng all branching ratios. This is achieved 
when the second argument of SetDecayChannels is set to true (it is false by 
default). 

In addition to the list of decay channels, a summary list of TTMDecay objects 
is generated in which each daughter appears only once, together with its total 
decay fraction. This summary list is automatically generated from the decay 
channel list when the SetDecayChannels function is called. 

An existing TList can be set as the decay channel list of the particle, using 
the SetDecayChannels function. This function calls UpdateDecaySummary, 
thereby automatically ensuring consistency between the decay channel and 
decay summary lists. 

The function SetDecayChannelEf f iciency sets the reconstruction efficiency 
of the specified decay channel to the specified percentage. Again, a consistent 
decay summary list is generated. 

Access to the TTMDecayChannel objects in the decay channel list is achieved 
through the GetDecayChannel method. If the extracted decay channel is sub- 
sequently altered, UpdateDecaySummary must be called to ensure consistency 
of the summary list. 



3.3 The TTMParticleSet Class 



The thermalised fireballs considered in statistical-thermal models typically 
contain approximately 350 different hadron and hadronic resonance species. 
To facihtate fast retrieval of particle properties, the TTMParticle objects of 
all constituents are stored in a hash table in a TTMParticleSet object. Other 
data members of this TTMParticleSet class include the filename used to in- 
stantiate the object and the number of particle species. Access to the entries 
in the hash table is through the particle Monte Carlo ID's. 
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3.3.1 Instantiating a TTMParticleSet Object 

In addition to the default constructor, the following constructors exist: 

TTMParticleSet *set = new TTMParticleSet (char *f iie) ; 
TTMParticleSet *set = new TTMParticleSet (TDatabasePDG *pdg) ; 

The first constructor instantiates a TTMParticleSet object and inputs the 
particle properties contained in the specified text file. As an example of such 
a file, /$THERMUS/particles/PartList_PPB2002.txt contains a list of all 
mesons (up to the -^'4(2045)) and baryons (up to the fl~) hsted in the July 
2002 Particle Physics Booklet [SB] (195 entries). Only particles need be in- 
cluded, since the anti-particle properties are directly related to those of the 
corresponding particle. The required file format is as follows: 

Delta(1600)0 32114 4 +1 1.60000 1 

0.35000 1.07454 (npiO) 

• stability flag (1 for stable, for unstable) 

• particle name 

• Monte Carlo particle ID (used for all referencing) 

• spin degeneracy 

• statistics (-1-1 for Fermi-Dirac, -1 for Bose-Einstein, for Boltzmann) 

• mass in GeV 

• strangeness 

• baryon number 

• charge 

• absolute strangeness content l^l 

• width in GeV 

• threshold in GeV 

• string recording the decay channel from which the threshold is calculated if 
the particle's width is non-zero 

All further particle properties have to be set with the relevant 'setters' (e.g. 
the charm, absolute charm content and hard-sphere radius). By default, all 
properties not listed in the particle list file are assumed to be zero. 

Figure E] shows the distribution of resonances (both particle and anti-particle) 
derived from /$THERMUS/particles/PartList_PPB2002 .txt. As collider en- 
ergies increase, so does the need to include also the higher mass resonances. 
Although the TTMParticle class allows for the properties of charmed parti- 
cles, these particles are not included in the default THERMUS particle list. 
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(A 160 




Mass (GeV) 

Fig. 3. The mass distribution of tlie resonances included in PartList_PPB2002.txt. 

If required, these particles have to be input by the user. The same applies to 
the hadrons composed of b and t quarks. 

It is also possible to use a TDatabasePDG object to instantiate a particle se 
TDatabasePDG objects also read in particle information from text files. The 
default file is /$ROOTSYS/etc/pdg_table .txt and is based on the parameters 
used in PYTHIA [57]. 

The constructor TTMParticleSet (TDatabasePDG *pdg) extracts only those 
particles in the specified TDatabasePDG object in particle classes 'Meson', 
'CharmedMeson', 'HiddenCharmMeson', 'B-Meson', 'Baryon', 'CharmedBaryon' 
and 'B-Baryon', as specified in /$ROOTSYS/etc/pdg_table.txt, and includes 
them in the hadron set. Anti-particles must be included in the TDatabasePDG 
object, as they are not automatically generated in this constructor of the 
TTMParticleSet class. 

The default file read into the TDatabasePDG object, however, is incomplete; the 
charm, degeneracy, threshold, strangeness, \S\, beauty and topness of the par- 
ticle are not included. Although the TDatabasePDG: : ReadPDGTable function 
and default file allow for isospin, I3, spin, flavor and tracking code to be entered 
too, the default file does not contain these values. Furthermore, all particles 



^ In order to have access to TDatabasePDG and related classes, one must first load 
/$ROOTSYS/lib/libEG . so 
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are made stable by default. Therefore, at present, using the TDatabasePDG 
class to instantiate a TTMParticleSet class should be avoided, at least until 
pdg_table.txt is improved. 



3.3.2 Inputting Decays 

Once a particle set has been defined, the decays to the stable particles in 
the set can be determined. After instantiating a TTMParticleSet object and 
settling on its stable constituents (the list of stable particles can be modified 
by adjusting the stability flags of the TTMParticle objects included in the 
TTMParticleSet object), decays can be input using the InputDecays method. 
Running this function populates the decay lists of all unstable particles in the 
set, using the decay files listed in the directory specified as the first argument. 
If a file is not found, then the corresponding particle is set to stable. For each 
typically unstable particle in /$THERMUS/particles/PartList_PPB2002 . txt, 
there exists a file in /$THERMUS/particles listing its decays. The filename 
is derived from the particle's name (e.g. Delta(1600) 0_decay.txt for the 
A(1600)°). There are presently 195 such files, with entries based on the Par- 
ticle Physics Booklet of July 2002 [56]. The decays of the corresponding 
anti-particles are automatically generated, while a private recursive function, 
GenerateBRatios, is invoked to ensure that only stable particles feature in 
the decay summary lists. The second argument of InputDecays, when set to 
true, scales the branching ratios so that their sum is 100%. As an example, 
consider the following (again only part of the listing is shown): 



root [ ] TTMParticleSet set("$THERMUS/particles/PartList_PPB2002.txt") 

root [ ] set . InputDecays("$THERMUS/particles/" , true) 

root [ ] TTMParticle *part = set . GetParticle (32114) 

root [ ] part->List() 



********* LISTING FOR PARTICLE Delta(1600)0 ********* 



UNSTABLE 



Decay Channels : 
BRatio: 0.108558 



Daughters : 
Daughters : 
Daughters : 
Daughters : 



2112 



111 



BRatio: 0.0542326 
BRatio: 0.272837 



2212 



2214 



-211 



-211 



BRatio: 0.0341395 



2114 



111 
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BRatio: 


0. 


.204651 


Daughters : 


1114 


211 


BRatio: 


0. 


. 0774884 


Daughters : 


2112 


113 


BRatio: 


0, 


.0387907 


Daughters : 


2212 


-213 


BRatio: 


0, 


. 139535 


Daughters : 


12112 


111 


BRatio: 


0, 


.0697674 


Daughters : 


12212 


-211 



Summary of Decays : 



2112 60.6774% 

111 62.5704% 

2212 39.3226% 

-211 83.9999% 

211 44.6773% 



For particle sets based on TDatabasePDG objects, decay lists should be pop- 
ulated through the function InputDecays (TDatabasePDG *). This function, 
however, does not automatically generate the anti-particle decays from those 
of the particle. Instead, the anti-particle decay list is used. Since the decay 
list may include electromagnetic and weak decays to particles other than the 
hadrons stored in the TTMParticleSet object, each channel is first checked 
to ensure that it contains only particles listed in the set. If not, the channel 
is excluded from the hadron's decay list used by THERMUS. As mentioned 
earlier, care should be taken when using TDatabasePDG objects based on the 
default file, as it is incomplete. 

An extremely useful function is ListParents (Int_t id) , which lists all of the 
parents of the particle with the specified Monte Carlo ID. This function uses 
GetParents(TList *parents, Int_t id), which populates the list passed 
with the decays to particle id. Note that these parents are not necessarily 
'direct parents'; the decays may involve unstable intermediates. 



3.3.3 Customising the Set 

The AddParticle and RemoveParticle functions allow customisation of par- 
ticle sets. Particle and anti-particle are treated symmetrically in the case of 
the former; if a particle is added, then its corresponding anti-particle is also 
added. This is not the case for the RemoveParticle function, however, where 
particle and anti-particle have to be removed separately. 

Mass-cuts can be performed using MassCut (Double_t x) to exclude all hadrons 
with masses greater than the argument (expressed in GeV) . Decays then have 
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to be re-inserted, to remove the influence of the newly-excluded hadrons from 
the decay lists. 

The function SetDecayEf f iciency allows the reconstruction efficiency of the 
decays from a specified parent to the specified daughter to be set. Changes 
arc reflected only in the decay summary list of the parent (i.e. not the decay 
channel list). Note that running UpdateDecaySummary or GenerateBRatios 
will remove any such changes, by creating again a summary list consistent 
with the channel list. 

In addition to these operations, users can input their own particle sets by 
compiling their own particle lists and decay files. 



3.4 The TTMParameter Class 

This class groups all relevant information for parameters in the statistical- 
thermal model. Data members include: 



f Name 


- the parameter name. 


fValue 


- the parameter value. 


f Error 


- the parameter error. 


fFlag 


- a flag signalling the type of parameter (constrain, flt. 




flxed, or uninitialised). 


f Status 


- a string reflecting the intended treatment or action taken 



In addition to these data members, the following, relevant to flt-type param- 
eters, are also included: 

f Start - the starting value in a fit, 

f Min - the lower bound of the fit-range, 

f Max - the upper bound of the fit-range, 

f Step - the step-size. 

The constructor and SetParameterdString name, Double_t value, Double_t 
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error) function set the parameter to fixed-type, by default. The parameter- 
type can be modified using the Constrain, Fit or Fix methods. 

3.5 The TTMParameterSet Class 

The TTMParameterSet class is the base class for all thermal parameter set 
classes. Each derived class contains its own TTMParameter array, with size 

determined by the requirements of the ensemble. The base class contains a 
pointer to the first element of this array. In addition, it stores the constraint 
information. 

All derived classes contain the function GetRadius. In this way, TTMParameterSet 
is able to define a function, GetVolume, which returns the volume required to 
convert densities into total fireball quantities. 

TTMParameterSetBSQ, TTMParameterSetBQ and TTMParameterSetCanBSQ are 

the derived classes. 



3.5.1 TTMParameterSetBSQ 

This derived class, applicable to the grand-canonical ensemble, contains the 
parameters: 

T hb I^s fJ^Q l^c Is Ic R, 

where R is the fireball radius, assuming a spherical fireball (i.e. V = A/SttR^). 
In addition, the B/2Q ratio and charm and strangeness density of the system 
are stored here. In the constructor, all errors are defaulted to zero, as is R, 
jic, S/V, C/V and B/2Q, while 70 is defaulted to unity. 

Each parameter has a 'getter' (e.g. GetTPar), which returns a pointer to 
the requested TTMParameter object. In this class, ns and hq can be set to 
constrain-type using ConstrainMuS and ConstrainMuQ, where the arguments 
are the required strangeness density and B/2Q ratio, respectively. No such 
function exists for /ic, since constraining functions are not yet coded for the 
charm density. Each parameter of this class can be set to fit-type, using func- 
tions such as FitT (where the fit parameters have reasonable default values), 
or fixed-type, using functions such as FixMuB. 
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3.5.2 TTMParameterSetBQ 



This derived class, applicable to the strangeness-canonical ensemble (strangeness 
exactly conserved and B and Q treated grand-canonically ) , has the parame- 
ters: 

T Hb I^q 7s Rc R, 

where Rc is the canonical or correlation radius; the radius inside which strangeness 
is exactly conserved. The fireball radius R, on the other hand, is used to con- 
vert densities into total fireball quantities. In addition, the required B/2Q 
ratio is also stored, as well as the strangeness required inside the correlation 
volume (which must be an integer). 

In addition to the same 'getters' and 'setters' as the previous derived class, it is 
possible to set /iq to constrain- type by specifying the B/2Q ratio in the argu- 
ment of ConstrainMuQ. The strangeness required inside the canonical volume 
is set through the SetS method. This value is defaulted to zero. The function 
ConserveSGlobally fixes the canonical radius, Rc, to the fireball radius, R. 
As in the case of the TTMParameterSetBSQ class, there also exist functions to 
set each parameter to fit or fixed-type. 



3.5.3 TTMParameterSetCanBSQ 

This set, applicable to the canonical ensemble with exact conservation of B, 
S and Q, contains the parameters: 

T B S Q -fs R. 

Since all conservation is exact, there arc no chemical potentials to satisfy 
constraints. Again, the same 'getters', 'setters' and functions to set each pa- 
rameter to fit or fixed- type exist, as in the case of the previously discussed 
TTMParameterSet derived classes. 



3.5.4 Example 

As an example, let us define a TTMParameterSetBQ object. By default, all 
parameters are initially of fixed-type. Suppose we wish to fit T and /ib, and 
use iiQ to constrain the B/2Q ratio in the model to that in Pb-|-Pb coUisions: 
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root [] TTMParameterSetBQ parBQ(0. 160,0.2,-0.01,0.8,6. ,6. ) 

root [ ] parBQ.FitKO. 160) 

root [ ] parBQ.FitMuB(0.2) 

root [ ] parBQ.ConstrainMuQC 1.2683) 

root [ ] parBQ.ListO 

***************************** Thermal Parameters **************************** 

Strangeness inside Canonical Volume = 

T = 0.16 (to be FITTED) 

start: 0.16 
range: 0.05 ~ 0.18 
step: 0.001 

(to be FITTED) 
start: 0.2 
range : — 0.5 
step: 0.001 

(to be CONSTRAINED) 

B/2q: 1.2683 

(FIXED) 

(FIXED) 

(FIXED) 



muB 



muQ 

gammas 
Can. radius 
radius 



0.2 



-0.01 

0.8 
6 
6 



Parameters unconstrained 



Note the default parameters for the T and jiB fits. Obviously, no constraining 
or fitting can take place yet; we have simply signalled our intent to take these 
actions at some later stage. 



3.6 The TTMThermalPartzcle Class 



By combining a TTMParticle and TTMParameterSet object, a thermal particle 
can be created. The TTMThermalParticle class is the base class from which 
thermal particle classes relevant to the three currently implemented thermal 
model formalisms, TTMThermalParticleBSQ, TTMThermalParticleBQ and 
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TTMThermalParticleCanBSQ, arc derived. Since no particle set is specified, the 
total fireball properties cannot be determined. Thus, in the grand-canonical 
approach, the constraints cannot yet be imposed to determine the values of 
the chemical potentials of constrain-type, while, in the strangeness-canonical 
and canonical formalisms, the canonical correction factors cannot yet be calcu- 
lated. Instead, at this stage, the chemical potentials and/or correction factors 
must be specified. 

Use is made of the fact that, in the Boltzmann approximation, E/V , N/V 
and P, in the canonical and strangeness-canonical ensembles, are simply the 
grand-canonical values, with the chemical potential(s) corresponding to the 
canonically-treated quantum number(s) set to zero, multiplied by a particle- 
specific correction factor. This allows the functions for calculating E/V, N/V 
and P in the Boltzmann approximation to be included in the base class, which 
then also contains the correction factor as a data member (by definition, this 
correction factor is 1 in the grand-canonical ensemble). 

Both functions including and excluding resonance width, F, are coded (e.g. 
DensityBoltzmannNoWidth and EnergyBoltzmannWidth). When width is in- 
cluded, a Breit-Wigner distribution is integrated over between the limits 
[max(m - 2T, mthreshow), + 21]. 



3.6.1 TTMThermalParticleBSQ 

This class is relevant to the grand-canonical treatment of B, S and Q. In ad- 
dition to the functions for calculating E/V, N/V and P in the Boltzmann ap- 
proximation, defined in the base class, functions implementing quantum statis- 
tics for these quantities exist in this derived class (e.g. EnergyQStatNoWidth 
and PressureQStatWidth). Additional member functions of this class calcu- 
late the entropy using either Boltzmann or quantum statistics, with or without 
width. 

In the functions calculating the thermal quantities assuming quantum statis- 
tics, it is first checked that the integrals converge for the bosons (i.e. there is no 
Bose- Einstein condensation). The check is performed by the ParametersAllowed 
method. A warning is issued if there are problems and zero is returned. 

This class also accommodates charm, since the associated parameter set in- 
cludes lie and 7c, while the associated particle may have non-zero charm. 
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3.6.2 TTMThermalParticleBQ 



This class is relevant to the strangeness-canonical ensemble. At present, this 
class is only applied in the Boltzmann approximation. Under this assumption, 
N/V, E/V and P are given by the grand-canonical result, with /is set to zero, 
up to a multiplicative correction factor. Since the total entropy does not split 
into the sum of particle entropies, no entropy calculation is made in this class. 



3.6.3 TTMThermalParticleCanBSQ 

This class is relevant to the fully canonical treatment oi B, S and Q. At 
present, as in the case of TTMThermalParticleBQ, this class is only apphed 
in the Boltzmann approximation. Also, since the total entropy again does not 
split into the sum of particle entropies, no entropy calculation is made here. 



3.6.4 Example 

Let us construct a thermal particle, within the strangeness-canonical ensemble, 
from the A(1600)° and the parameter set previously defined. Since this particle 
has zero strangeness, a correction factor of 1 is passed as the third argument 
of the constructor: 

root [ ] TTMThermalParticleBQ thenii_delta(part ,&parBQ , 1 . ) 

root [ ] therm_delta.DensityBoltzmaimNoWidth() 

(Double_t) 8.150726717100899136-04 

root [ ] theriii_delta.EnergyBoltzmannWidth() 

(Double_t) 2 . 29185316377137748e-03 

3.7 The TTMThermalModel Class 

Once a parameter and particle set have been specified, these can be com- 
bined into a thermal model. TTMThermalModel is the base class from which 
the TTMThermalModelBSQ, TTMThermalModelBQ and TTMThermalModelCanBSQ 
classes are derived. A string descriptor is included as a data member of the 
base class to identify the type of model. This is used, for example, to handle 
the fact that the number of parameters in the associated parameter sets is 
different, depending on the model type. 

All derived classes define functions to calculate the primordial particle, energy 
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and entropy densities, as well as the pressure. These thermal quantities are 
stored in a hash table of TTMDensObj objects. Again, access is through the 
particle ID's. In addition to the individual particles' thermal quantities, the 
total primordial fireball strangeness, baryon, charge, charm, energy, entropy, 
and particle densities, pressure, and Wroblewski factor (see Section [3. 7. Ill) are 
included as data members. 

At this level, the constraints on any chemical potentials of constrain-type can 
be imposed, and the correction factors in canonical treatments can be deter- 
mined. Also, as soon as the primordial particle densities are known, the decay 
contributions can be calculated. 



3. 7. 1 Calculating Particle Densities 

Running GenerateParticleDens clears the current entries in the density hash 
table of the TTMThermalModel object, automatically constrains the chem- 
ical potentials (where applicable), calculates the canonical correction fac- 
tors (where applicable), and then populates the density hash table with a 
TTMDensObj object for each particle in the associated set. The decay contribu- 
tions to each stable particle are also calculated, so that the density hash table 
contains both primordial and decay particle density contributions, provided of 
course that the decays have been entered in the associated TTMParticleSet 
object. In addition, the Wroblewski factor and total strangeness, baryon, 
charge, charm and particle densities in the fireball are calculated. 

Note: The summary decay lists of the associated TTMParticleSet object are 
used to calculate the decay contributions. Hence, only stable particles have 
decay contributions reflected in the hash table. Unstable particles that are 
themselves fed by higher-lying resonances, do not receive a decay contribution. 

Each derived class contains the private function PrimPartDens, which calcu- 
lates only the primordial particle densities and, hence, the canonical correction 
factors, where applicable. In the case of the grand-canonical and strangeness- 
canonical ensembles, this function calculates the densities without automat- 
ically constraining the chemical potentials of constrain-type first. The con- 
straining is handled by GenerateParticleDens, which calls external friend 
functions, which, in turn, call PrimPartDens. In the purely canonical ensem- 
ble, GenerateParticleDens simply calls PrimPartDens. In this way, there is 
uniformity between the derived classes. Since there is no constraining to be 
done, there is no real need for a separate function in the canonical case. 
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3.1.2 Calculating Energy and Entropy Densities and Pressure 

GenerateEnergyDens, GenerateEntropyDens and GeneratePressure iterate 
through the existing density hash table and calculate and insert, respectively, 
the primordial energy density, entropy density and pressure of each particle 
in the set. In addition, they calculate the total primordial energy density, 
entropy density and pressure in the fireball, respectively. These functions re- 
quire that the density hash table already be in existence. In other words, 
GenerateParticleDens must already have been run. If the parameters have 
subsequently changed, then this function must be run yet again to recalculate 
the correction factors or re-constrain the parameters, as required. 



3.7.3 Bose-Einstein Condensation 

When quantum statistics are taken into account (e.g. in TTMThermalModelBSQ 
or for the non-strange particles in TTMThermalModelBQ), certain choices of pa- 
rameters lead to diverging integrals for the bosons (Bose-Einstein condensa- 
tion). In these classes, a check, based on TTMThermalParticleBSQ : : Parameters- 
Allowed, is included to ensure that the parameters do not lead to problems. 
Including also the possibility of incomplete strangeness and/or charm satu- 
ration (i.e. 75 7^ 1 and/or jc !)> Bose-Einstein condensation is avoided, 
provided that, 

g(mi-Mi)/T ^ ^^5i|^JC,|^ (38) 



for each boson. If this condition is failed to be met for any of the bosons in 
the set, a warning is issued and the densities are not calculated. 



3.7.4 Accessing the Thermal Densities 

The entries in the density hash table are accessed using the particle Monte 
Carlo id's. The function GetDensities(Int_t ID) returns the TTMDensObj 
object containing the thermal quantities of the particle with the specified ID. 
The primordial particle, energy, and entropy densities, pressure, and decay 
density are extracted from this object using the GetPrimDensity, GetPrimEnergy, 
GetPrimEntropy, GetPrimPressure, and GetDecayDensity functions of the 
TTMDensObj class, respectively. The sum of the primordial and decay particle 
densities is returned by TTMDensObj : :GetFinalDensity. TTMDensObj : :List 
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outputs to screen all thermal densities stored in a 
TTMDensObj object. 

ListStableDensities lists the densities (primordial and decay contributions) 
of all those particles considered stable in the particle set associated with the 
model. Access to the total fireball densities is through separate 'getters' de- 
fined in the TTMThermalModel base class (e.g. GetStrange, GetBaryon etc.). 



3. 7. 5 Further Functions 

GenerateDecayPartDens and GenerateDecayPartDens(Int_t id) (both de- 
fined in the base class) calculate decay contributions to stable particles. The 
former iterates through the density hash table and calculates the decay contri- 
butions to all those particles considered stable in the set. The latter calculates 
just the contribution to the stable particle with the specified ID. In both cases, 
the primordial densities must be calculated first. In fact, GenerateParticleDens 
automatically calls GenerateDecayPartDens, so that this function does not 
have to be run separately under ordinary circumstances. However, if one is 
interested in investigating the effect of decays, while keeping the parameters 
(and hence the primordial densities) fixed, then running these functions is best 
(the hash table will not be repeatedly cleared and repopulated with the same 
primordial densities). 

ListDecayContributions(Int_t d_id) lists the contributions (in percent- 
age and absolute terms) of decays to the daughter with the specified ID. The 
primordial and decay densities must already appear in the density hash ta- 
ble (i.e. run GenerateParticleDens first). ListDecayContribution(Int_t 
p_id, Int_t d_id) lists the contribution of the decay from the specified par- 
ent (with ID p_id) to the specified daughter (with ID d_id). The percentages 
listed by each of these functions are those of the individual decays to the total 
decay density. 

Next we consider the specific features of the derived TTMThermalModel classes. 

3. 7. 6 TTMTherma IMode IBSQ 

In the grand-canonical ensemble, quantum statistics can be employed and, 
hence, there is a fiag specifying whether to use Fermi-Dirac and Bose-Einstein 
statistics or Boltzmann statistics. The constructor, by default, includes both 
the effect of quantum statistics and resonance width. The flags controUing their 
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inclusion are set using the SetQStats and SetWidth functions, respectively. 
The functions that calculate the particle, energy, and entropy densities, and 
pressure then use the corresponding functions in the TTMThermalParticleBSQ 
class to calculate these quantities in the required way. The statistics data mem- 
ber (f Stat) of each TTMParticle included in the associated set can be used 
to fine-tune the inclusion of quantum statistics; with the quantum statistics 
fiag switched on, Boltzmann statistics are still used for those particles with 
fStat=0. 

In this ensemble, at this stage, both /i^ and fig can be constrained (either 
separately or simultaneously). In order to accomplish this, the /ig and/or 
fiQ parameters in the associated TTMParameterSetBSQ object must be set to 
constrain-type. 

It is also possible to constrain fis by the primordial ratio E/N (the average 
energy per hadron), nj^ + ni (the total primordial baryon plus anti-baryon den- 
sity), or s/T^ (the primordial, temperature- normalised entropy density). This 
is accomplished by the ConstrainEoverN, ConstrainTotalBaryonDensity 
and ConstrainSoverTS methods, respectively. Running these functions will 
adjust /is such that E /N , Ub + ni or s/T^, respectively, has the required value, 
regardless of the parameter type of /x^. In addition, the percolation model \5E\ 
can be imposed to constrain /^^ using ConstrainPercolation. 

This class also accommodates charm, since the associated parameter set in- 
cludes fic and 7c, while the associated particle set may contain charmed 
particles. However, no constraining functions have yet been written for the 
charm content within this ensemble. 

Within the grand-canonical ensemble, it is possible to include excluded vol- 
ume effects. Their inclusion is controlled by the f ExclVolCorrection flag, 
false by default, which is set through the SetExcludedVolume function. When 
included, these corrections are calculated on calling GenerateParticleDens, 
based on the hard-sphere radii stored in the TTMParticle objects of the asso- 
ciated particle set. 



3.7.7 TTMThermalModelBQ 

This class contains the following additional data members: 
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X ~l «. "74- ^ 4- 

ilnZtot 


- log of the total partition function, 


I InZO 


- log of the non-strange component of the partition function, 


f ExactMuS 


- equivalent strangeness chemical potential. 


f CorrPl 


- canonical correction for S = +1 particles. 


f CorrP2 


- canonical correction for S = +2 particles. 


fCorrPS 


- canonical correction for S = +3 particles. 


fCorrMl 


- canonical correction for S = -1 particles, 


fCorrM2 


- canonical correction for S = -2 particles. 


fCorrMS 


- canonical correction for S = -3 particles. 



Although this ensemble is only applied in the Boltzmann approximation for 
S ^ hadrons, it is possible to apply quantum statistics to the S* = hadrons. 
This is achieved through the SetNonStrangeQStats function. By default, 
quantum statistics is included for the non-strange hadrons by the construc- 
tors. Resonance width can be included for all hadrons, and is achieved through 
the SetWidth function. The constructors, by default, apply resonance width. 
The functions that calculate the particle, energy, and entropy densities, and 
pressure then use the corresponding functions in the TTMThermalParticle 
classes to calculate these quantities in the required way. 

GenerateParticleDens populates the density hash table with particle den- 
sities, including the canonical correction factors, which are also stored in the 
appropriate data members. The equivalent strangeness chemical potential is 
calculated from the canonical correction factor for 5" = -|-1 particles. In the 
limit of large VT"^, this approaches the value of /is in the equivalent grand- 
canonical treatment. 



Running GenerateEntropyDens populates each TTMDensObj object in the 
hash table with only that part of the total entropy that can be unambiguously 
attributed to that particular particle. There is a term in the total entropy that 
cannot be split; this is added to the total entropy at the end, but not included 
in the individual entropies (i.e. summing up the entropy contributions of each 
particle will not give the total entropy). 

At this stage, in this formalism, fig can be constrained (this is automati- 
cally realised if this parameter is set to constrain-type), while the correlation 
radius (i?c) can be set to the fireball radius (i?) by applying the function 
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ConserveSGlobally to the associated TTMParameterSetBQ object. 

In exactly the same way as in the grand-canonical ensemble case, iib can be 
constrained in this ensemble by the primordial ratio E/N (the average energy 
per hadron), rib + ni (the total primordial baryon plus anti-baryon density), 
or s/T^ (the primordial, temperature-normalised entropy density), as well as 
by the percolation model. 

3. 7. 8 TTMTherma IModelCanBSQ 

This class contains, amongst others, the following data members: 
f InZtot - log of the total canonical partition function. 



-fMiiC -FMiiQ -fMnn 


- equivalent 


chemical potentials. 


^^^^^^ T 1^ 

iL>orrpip 


- correcxion 


for 


71 -iiKe parxicies. 


f Corrpim 


- correction 


for 


7r~-like particles. 


f Corrkm 


- correction 


for 


K~-\ike particles. 


f Corrkp 


- correction 


for 


i^'^-like particles. 


fCorrkO 


- correction 


for 


_ft'*'-like particles. 


fCorrakO 


- correction 


for 


K^-\ike particles. 


f Corrproton 


- correction 


for 


p-like particles. 


f Corraproton 


- correction 


for 


p-like particles. 


f Corrneutron 


- correction 


for 


n-like particles. 


f Corr aneut r on 


- correction 


for 


n-like particles. 


f Corrlambda 


- correction 


for 


A-like particles. 


f Corralambda 


- correction 


for 


A-hke particles. 


f Corrsigmap 


- correction 


for 


E+-hke particles. 


f Corrasigmap 


- correction 


for 


S~-like particles. 


f Corrsigmam 


- correction 


for 


E~-like particles. 


f Corrasigmam 


- correction 


for 


E+-hke particles. 


f Corrdeltam 


- correction 


for 


A~-like particles. 
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f Corradeltam 


- correction 


for 


A"'"-like particles, 


f Corrdeltapp 


- correction 


for 


A^~''-like particles, 


f Corradeltapp 


- correction 


for 


A -like particles. 


f Corrksim 


- correction 


for 


S~-like particles. 




- correcxion 


for 


^ -iiKe parxicies. 


fCorrksiO 


- correction 


for 


S°-like particles. 


f CorraksiO 


- correction 


for 


S*^-like particles, 


fCorr omega 


- correction 


for 


f]~-like particles. 


f Corraomega 


- correction 


for 


Q+-like particles. 



Since this ensemble is only applied in the Boltzmann approximation, there is 
no flag for quantum statistics. However, resonance width can be included. This 
is achieved through the SetWidth function. The constructor, by default, ap- 
plies resonance width. The functions that calculate the particle, energy, and 
entropy densities, and pressure then use the corresponding functions in the 
TTMThermalParticle classes to calculate these quantities in the required way. 

GenerateParticleDens calls PrimPartDens, which calculates the particle den- 
sities, including the canonical correction factors, which are then also stored 
in the relevant data members accessible through the GetCorrFactor method. 
The integrands featuring in the evaluation of the partition function and correc- 
tion factors can be viewed after calling PopulateZHistograms. This function 
populates the array passed as argument with histograms showing these in- 
tegrands as a function of the integration variables (ps and Since these 
histograms are created off of the heap, they must be cleaned up afterwards. 

GenerateEntropyDens acts in exactly the same way as in the strangeness- 
canonical ensemble case. 



3.7.9 Example 

As an example, we consider the strangeness-canonical ensemble, based on the 
particle set and strangeness-canonical parameter set previously defined. After 
instantiating the object, we populate the hash table with primordial and decay 
particle densities: 
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root [ ] TTMThermalModelBQ modBQ(&set ,&parBQ) 
root [ ] modBQ .GenerateParticleDensO 
root [ ] parBQ.ListO 



***************************** Thermal Parameters **************************** 



Strangeness inside Canonical Volume = 



muB 



muQ 

gammas 
Can. radius 
radius 



0.16 



0.2 



-0.00636409 

0.8 
6 



(to be FITTED) 
start: 0.16 
range: 0.05 ~ 0.18 
step: 0.001 

(to be FITTED) 
start: 0.2 
range : — 0.5 
step: 0.001 

(♦CONSTRAINED*) 

B/2Q: 1.2683 

(FIXED) 

(FIXED) 

(FIXED) 



B/2Q Successfully Constrained 



One notices that the constraint on /iq is now automatically imposed. 



The energy and entropy densities and pressure can be calculated once GenerateParticleDens 
has been run: 

root [ ] modBQ. GenerateEnergyDensO 
root [ ] modBQ. GenerateEntropyDensO 
root [ ] modBQ. GeneratePressureO 

Now, suppose that we are interested in the thermal densities of the A(1600)'' 
and TT"*": 



root [ ] TTMDensObj *delta_dens = modBQ. GetDensities (32114) 
root [ ] delta_dens->List() 
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**** Densities for Particle 32114 **** 
n_prim = 0.00138306 
n_decay = 
e_prim = 0.0022912 
s_prim = 0.0139745 
p_prim = 0.000221328 

root [ ] TTMDensDbj *piplus_dens = modBQ . GetDensities (211) 
root [ ] piplus_dens->List 

**** Densities for Particle 211 **** 

n_prim = 0.0488139 

n_decay = 0.119683 

e_prim = 0.0247039 

s_prim = 0.20276 

p_prim = 0.00742708 



One notices that the tt"*" has a decay density contribution, while the A(1600)° 
does not. This is because, unhke the A(1600)°, the tt"*" was considered stable. 



3.7.10 Imposing of Constraints 

The 'Numerical Recipes in C [59] function applying the Broyden globally con- 
vergent secant method of solving nonlinear systems of equations is employed 
by THERMITS to constrain parameters. The input to the Broyden method 
is a vector of functions for which roots are sought. Typically, in the thermal 
model, solutions to the following equations are required (either separately or 
simultaneously) : 



CT-> /T r \ model / T-> \ collidinq system 

-(pi =0, 
'^Q J primordial ) 

nmodel ncolliding system n 

'-^primordial '-^ ^i 



primordial 

lod 
'in 

model / TP \ required 



E/V\ [E 



N/V) 



primordial 



N 



0. 



Although, as written, these equations are correct, the quantities -B/2Q, S 
and E/N are typically of different orders of magnitude. Since the Broyden 
method in 'Numerical Recipes in C defines just one tolerance level for function 
convergence (TOLF), it is important to 'normalise' each equation: 
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B/V 
2Q/V 



) 



primordial 



model 




This is the most democratic way of treating the constraints. However, this 
method obviously fails in the event of one of the denominators being zero. For 
the equations considered above, this is only likely in the case of the strangeness 
constraint, where the initial strangeness content is typically zero. In this case, 
where the strangeness carried by the positively strange particles 5*+ is balanced 
by the strangeness carried by the negatively strange particles S*., we write as 
our function to be satisfied, 

( ^ l-frynodel / { \ Q \model i-tr i i o model f) 

\yO I V Jpj-imordial I I I + I pr imordiai / ^ ' l^^lprimordiall 'J ^" 

In this way, the constraints can be satisfied to equal relative degrees, and 
equally well fractionally at each point in the parameter space. In addition 
to the constraints listed above, THERMUS also allows for the constraining 
of the total baryon plus anti-baryon density and the temperature-normalised 
entropy density, s/T^, as well as the imposing of the percolation model. 



3.7.11 Calculation of the Wrdblewski Factor 
The Wroblewski factor [60] is defined as. 



where < uu > + < dd > is the sum of newly-produced uu and dd pairs, while 
all ss pairs are newly-produced if S = in the initial state. 

In THERMUS, is calculated in the following way: 

• Using the primordial particle densities and the strangeness content of each 
particle listed in the particle hash table, the s + s and u + d + u + d densities 
are determined. 

• Assuming S = 0, #s = and so the density of newly-produced ss pairs 



A5 



2 < ss > 



< uu > + < dd >' 
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is simply (s + s)/2. 

• From baryon number conservation, the net baryon content in the system, 
hb-i originates from the initial state. Thus, 3 x must correspond to the 
density oiu + d quarks brought in by the coUiding nuclei. This is subtracted 
from the total u -\- d -\- u -\- d density to yield the density of newly-produced 
non-strange light quarks. 

• Since #s = #s and, amongst newly-produced non-strange light quarks, 
u + d = u + d, further assuming that = implies that u = u = d = d. 
This allows the density of uu and dd pairs to be easily determined. 

3.8 Thermal Fit Preliminaries 

Often a single experiment releases yields and ratios that contain different 
feed-down corrections. Each yield or ratio then has a different decay chain as- 
sociated with it. Since TTMThermalModel objects allow for just one associated 
particle set, they do not allow sufficient flexibility for performing thermal fits 
to experimental data. However, TTMThermalFit classes do feature such flex- 
ibility. Before we discuss these classes, let us look at the TTMYield object, 
which forms an essential part of the TTMThermalFit class. 



3.8.1 TTMYield 

Information relating to both yields and ratios of yields can be stored in 
TTMYield objects. These objects contain the following data members: 



f Name 


- the name of the yield or ratio. 


fIDl 


- the ID of the yield or numerator ID in the case of a ratio. 


f ID2 


- denominator ID in the case of a ratio (0 for a yield). 


fPit 


- true if the yield or ratio is to be included in a flt (else predicted). 


fSetl 


- particle set relevant to yield or numerator in case of ratio. 


fSet2 


- particle set relevant to denominator in case of ratio (0 for yield). 


f ExpValue 


- the experimental value. 


f ExpError 


- the experimental error. 


fModelValue 


- the model value. 


fModelError 


- the model error. 
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By default, TTMYield objects are set for inclusion in fits. The functions Fit 
and Predict control the fit-status of a TTMYield object. Particle sets (decay 
chains) are assigned using the SetPartSet method. 



The functions GetStdDev and GetQuadDev return the number of standard and 
quadratic deviations between model and experimental values, respectively, i.e.. 



respectively, while List outputs the contents of a TTMYield object to screen. 
Access to all remaining data members is through the relevant 'getters' and 
'setters'. 



3.9 The TTMThermalFit Class 

This is the base class from which the TTMThermalFitBSQ, TTMThermalFitBQ 
and TTMThermalFitCaiiBSQ classes are derived. Each TTMThermalFit object 
contains: 

• a particle set, the so-called base set, which contains all of the constituents 
of the hadron g well as the default decay chain to be used; 

• a parameter set; 

• a list of TTMYield objects containing yields and/or ratios of interest; 

• data members storing the total quadratic deviation; and 

• a TMinuit fit object. 

A string descriptor is also included in the base class to identify the type of 
model on which the fit is based. This is used, for example, to determine the 
number of parameters in the associated parameter sets. 

Each derived class defines a private function, GenerateThermalModel, which 

creates (off the heap) a thermal model object, based on the base particle set 
and parameter set of the TTMThermalFit object, with the specific quantum 
statistics/resonance width/excluded volume requirements, where applicable. 



(Model Value - Exp. Value) /Exp. Error, 



(39) 



and. 



(Model Value - Exp. Value) /Model Value, 



(40) 
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3. 9. 1 Populating and Customising the List of Yields of Interest 

The list of yields and/or ratios of interest can be input from file using the 
function InputExpYields, provided that the file has the following format: 



333 


Exp_A 


0.02 




0.01 






-211 


211 


Exp_B 


0, 


.990 


0, 


,100 


-211 


211 


Exp_C 


0, 


.960 


0, 


.177 


321 


-321 


Exp_C 


1, 


.152 


0, 


.239 



where the first line corresponds to a yield, and has format: 

Yield ID /t Descriptor string /t Exp. Value /t Exp. Error /n 

while the remaining lines correspond to ratios, and have format: 

Numerator ID /t Denominator ID /t Descriptor string /t Exp. Value 
/t Exp. Error/n 

A TTMYield object is created off the heap for each line in the file, with a 
name derived from the ID's and the descriptor. This name is determined by 
the private function GetName, which uses the base particle set to convert the 
particle ID's into particle names and appends the descriptor. In addition to 
all of the Monte Carlo particle ID's in the associated base particle set, the 
following THERMUS-defined identifiers are also allowed: 

• ID = 1: Np^ru 

• ID = 2: h-, 

• ID = 3: h+. 

A TTMYield object can also be added to the list using AddYield. Such yields 
should, however, have names that are consistent with those added by the 
InputExpYields method; the GetName function should be used to ensure this 
consistency. Only yields with unique names can be added to the list, since it 
is this name which allows retrieval of the TTMYield objects from the list. If a 
yield with the same name already exists in the list, a warning is issued. The 
inclusion of descriptors ensures that TTMYield objects can always be given 
unique names. 

RemoveYield(Int_t idl,Int_t id2,TString descr) removes from the list 
and deletes the yield with the name derived from the specified ID's and de- 
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scriptor by GetName. The GetYield(Int_t idl,Int_t id2,TString descr) 
method returns the required yield. 



3.9.2 Generating Model Values 

Values for each of the yields of interest hsted in a TTMThermalFit object are 
calculated by the function GenerateYields. This method uses the current 
parameter values and assigned particle sets to calculate these model values. 



GenerateYields firstly calculates the primordial particle densities of all con- 
stituents listed in the base particle set. This it does by creating the relevant 
TTMThermalModel object from the base particle set and the parameters, and 
then calling GenerateParticleDens. In this way, the density hash table of 
the newly-formed TTMThermalModel object is populated with primordial den- 
sities, as well as decay contributions, according to the base particle set (recall 
that GenerateParticleDens automatically calculates decay contributions in 
addition to primordial ones). GenerateYields then iterates through the list of 
TTMYield objects, calculating their specific decay contributions. New model 
values are then inserted into these TTMYield objects. In addition, the total 
and quadratic deviation are calculated, based solely on the TTMYield objects 
which are of fit-type. ListYields lists all TTMYield objects in the list. 



3.9.3 Performing a Fit 

The FitData(Int_t flag) method initiates a fit to all experimental yields 
or ratios in the TTMYield hst which are of fit-type. With f lag=0, a fit 
is performed, while f lag=l leads to a quadratic deviation fit. In both cases, 
fit_f unction is called. This function determines which parameters of the 
associated parameter set are to be fit, and performs the required fit using 
the ROOT TMinuit fit class. On completion, the list of TTMYield objects 
contains the model values, while the parameter set reflects the best-flt param- 
eters. Model values are calculated by the GenerateYields method. For each 
TTMYield object in the list, a model value is calculated- even those that have 
been chosen to be excluded from the actual fit. In this way, model predictions 
can be determined at the same time as a fit is performed. ListMinuitlnf o 
lists all information relating to the TMinuit object, following a fit. 
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3.9.4 TTMThermalFitBSQ, TTMThermalFitBQ and TTMThermalFitCanBSQ 



The constructor in each of these derived classes instantiates an object with 
the specified base particle set and parameter set and inputs the yields listed 
in the specified file in the TTMYield list. 

The specifics of the fit, i.e. the treatment of quantum statistics (in the grand- 
canonical ensemble and for the non-strange particles in the strangeness-canonical 
ensemble) , resonance width (in all three ensembles) and excluded volume cor- 
rections (in the grand-canonical ensemble), are handled through the SetQStats/ 
SetNonStrangeQStats, SetWidth and SetExclVol methods, respectively. By 
default, both resonance width and quantum statistics are included, while ex- 
cluded volume corrections are excluded, where applicable. 



3. 9. 5 Example 

As an example to conclude this section, consider a fit to fictitious particle 
ratios measured in Au-|-Au collisions at some energy. Wc will assume a grand- 
canonical ensemble, with the parameters T, and ns fitted, and /xq fixed to 
zero. In the grand-canonical ensemble, ratios are independent of the fireball 
radius (this is not true in the canonical ensemble). For this reason, there is no 
need to specify the treatment of the radius. Furthermore, we will ignore the 
effects of resonance width and quantum statistics. 

We begin by instantiating a particle set object, based on the particle list 
distributed with THERMUS. After inputting the particle decays (scaled to 
100%), a parameter set is defined: 

root [ ] TTMParticleSet set("$THERMUS/particles/PartList_PPB2002.txt") 
root [ ] set. IiiputDecays("$THERMUS/particles/", true) 
root [] TTMParameterSetBSQ par (0.160, 0.05,0. ,0. ,1.) 

Next, we change the parameters T, fis and fis to fit-type, supplying sensible 
starting values as the arguments to the appropriate functions. For all other 
properties of the fit (step size, fit range etc.), we accept the default values: 

root [ ] par.FitT(0. 160) 
root [ ] par. FitMuB (0.05) 
root [ ] par.FitMuS(0.) 

Next, we prepare a file ('ExpData.txt') containing the experimental data: 
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0. 
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0. 


.050 


-3122 


3122 


Exp. 


.B 


0. 


,734 


0. 


.210 


-3122 


3122 


Exp. 


.C 


0. 


,720 


0. 


.024 


-3312 


3312 


Exp. 


.C 


0, 


.878 


0, 


.054 


-3334 


3334 


Exp. 


.c 


1. 


,062 


0. 


.410 



As one can see, there are multiple occurrences of the same particle-anti- 
particle combination. This is why additional descriptors are required. In this 
case, the descriptors list the particular experiment responsible for the mea- 
surement. In other situations, the descriptors may describe whether feed-down 
corrections have been employed or some other relevant detail that, together 
with the id's, uniquely identifies each yield or ratio. 

We are now in a position to create a TTMThermalFitBSQ object based on 
the newly-instantiated parameter and particle sets and the data file. Since 
quantum statistics and resonance width are included by default, we have to 
explicitly turn these settings off: 

root [ ] TTMThermalFitBSQ f it (&set,&par, "ExpData.txt") 
root [ ] fit.SetQStats(kFALSE) 
root [ ] fit.SetWidth(kFALSE) 

Next, let us simply generate the model values corresponding to each of the 
TTMYield objects in the list, based on the current parameters. Part of the 
output of ListYields is shown here: 

root [ ] f it . GenerateYields () 
root [ ] fit. ListYields 



K+/anti-K+ Exp_B: 

FIT YIELD 

Experiment: 1.152 +- 0.239 
Model: 0.979833 +- 

Std.Dev.: -0.720365 Quad. Dev. : -0.175711 
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K+/anti-K+ Exp_D: 

FIT YIELD 

Experiment: 1.098 +- 0.111 
Model: 0.979833 +- 

Std.Dev.: -1.06457 Quad. Dev. : -0.120599 

K+/anti-K+ Exp_C: 

FIT YIELD 

Experiment: 1.108 +- 0.022 
Model: 0.979833 +- 

Std.Dev.: -5.82578 Quad. Dev. : -0.130805 

anti-p/p Exp_A: 

FIT YIELD 

Experiment: 0.65 +- 0.092 

Model: 0.535261 +- 

Std.Dev.: -1.24716 Quad. Dev. : -0.21436 



Finally, we perform a fit: 
root [ ] fit.FitData(O) 



COVARIANCE MATRIX CALCULATED SUCCESSFULLY 

FCN=3. 54326 FROM MIGRAD STATUS=CONVERGED 128 CALLS 129 TOTAL 

EDM=6.43919e-06 STRATEGY= 1 ERROR MATRIX ACCURATE 



EXT PARAMETER 
NO . NAME 

1 T 

2 muB 

3 muS 
EXTERNAL ERROR MATRIX 

4.593e-03 1.289e-03 



VALUE 

1.628786-01 
3.589088-02 
1.06828e-02 

NDIM= 25 
5.418e-04 



ERROR 

1.10211e-01 
1.913648-02 
8.13945e-03 
NPAR= 



1.2898-03 
5.418e-04 



3.6898-04 
1.548e-04 



1.5488-04 
6.6538-05 



STEP 

SIZE 
2.32021e-04 
1.685438-05 
1.80744e-05 



FIRST 
DERIVATIVE 
-7.988098-03 
4.057408-03 
1.25485e-01 



ERR DEF=1 
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PARAMETER CORRELATION COEFFICIENTS 

NO. GLOBAL 12 3 

1 0.99017 1.000 0.990 0.980 

2 0.99410 0.990 1.000 0.988 

3 0.98814 0.980 0.988 1.000 

FCN=3. 54326 FROM MIGRAD STATUS=CONVERGED 128 CALLS 129 TOTAL 

EDM=6.43919e-06 STRATEGY= 1 ERROR MATRIX ACCURATE 

PHYSICAL LIMITS 
ERROR NEGATIVE POSITIVE 

1.10211e-01 5.00000e-02 1.80000e-01 
1.91364e-02 O.OOOOOe+00 5.00000e-01 
8.13945e-03 O.OOOOOe+00 5.00000e-01 



EXT PARAMETER 
NO . NAME 

1 T 

2 muB 

3 muS 



VALUE 
1.62878e-01 
3.58908e-02 
1.06828e-02 



Once completed, the associated parameter set contains the best-fit values for 
the fit parameters: 

root [ ] par. List () 

***************************** Thermal Parameters **************************** 



0.162878 



0.110211 



(FITTED ! ) 
start: 0.16 
range : . 05 - 
step: 0.001 



0.18 



muB 



0.0358908 



0.0191364 



(FITTED ! ) 
start: 0.05 
range : — 0.5 
step: 0.001 



muS 



0.0106828 +- 0.00813945 



(FITTED ! ) 
start : 
range : — 0.5 
step: 0.001 



muQ 
gammas 
radius 

muC 
gammac 




1 



1 



(FIXED) 
(FIXED) 
(FIXED) 
(FIXED) 
(FIXED) 



Parameters unconstrained 
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4 Installation of THERMUS 



Having introduced the basic functionality of THERMUS in the previous sec- 
tion, we conclude by outlining the installation procedure. 

Since several functions in THERMUS use 'Numerical Recipes in C code [59] 
(which is under copyright), it is required that THERMUS users have their own 
copies of this software. Then, with ROOT [9j already installed, the following 
steps are to be followed to install THERMUS: 

• Download the THERMUS source; 

• Set an environment variable 'THERMUS' to point at the top-level directory 
containing the THERMUS code; 

• Copy the following 'Numerical Recipes in C [59] functions to $ (THERMUS) /nrc: 



broydn . 


c 


rsolv . c 


f dj ac . c 


fmin. c 


Insrch. 


c 


nrutil . c 


nrutil . 


h 


qrdcmp . c 


qrupdt . 


c 


rotate . c 


zbrent . 


c; 





• Use the makefiles in $ (THERMUS) /functions, $ (THERMUS) /nrc and $ (THERMUS) / 
to build the libFunctions . so, libNRCFunctions . so and libTHERMUS . so 
shared object files (run make all in each of these directories); 

• Finally, open a ROOT session, load the libraries and begin: 

root [ ] gSystem->Load(" ./lib/libFunctions.so") ; 
root [] gSystem->Load(" ./lib/libNRCFunctions.so") ; 
root [] gSystem->Load(" ./lib/libTHERMUS.so") ; 
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